Introduction

Ce programme permet de préparer l’ensemble des données statistiques et cartographiques utiles pour l’étude d’une agglomération urbaine à partir du fichier ménage du recensement de population 2019. Il charge de très gros fichiers contenus dans un dossier temporaire tmp et en extrait les informations utiles pour les placer dans un dossier data/nomagg . Il a besoin en entrée du code et du nom de la commune centre ainsi que du chemin vers le dossier de stockage des résultats

# Commune centre 
codectr<-  "35238"  # Code
namectr <- "Rennes"  # Nom

# Dossier de stockage
myrep <- "data/Rennes/"

# Choix du fichier ménage et de ses métadonnées
ficmen<-"tmp/FD_LOGEMTZC_2019.csv"
ficmen_meta <-"tmp/Varmod_LOGEMT_2019.csv"

Etape 1 : Données administratives

On commence par récupérer les informations adminstratives sur la zone d’étude

adm<-read.table("tmp/maps/table-appartenance-geo-communes-23.csv",
                sep=";", 
                header=T, fileEncoding = "UTF-8",
                quote = '"')
head(adm)

adm_ctr<-adm %>% filter(adm$CODGEO==codectr)
adm_ctr

my_EPCI <- adm_ctr$EPCI
my_AA <- adm_ctr$AAV2020

AA_adm <- adm %>% filter(AAV2020 == my_AA)
saveRDS(AA_adm, paste0(myrep,"AA_adm.RDS"))

EPCI_adm <- adm %>% filter(EPCI == my_EPCI)
saveRDS(EPCI_adm, paste0(myrep,"EPCI_adm.RDS"))

CTR_adm <- adm %>% filter(CODGEO == codectr)
saveRDS(CTR_adm, paste0(myrep,"CTR_adm.RDS"))

Etape 2 : Données géométriques

On va maintenant extraire les données géométriques relatives au contour des communes et des IRIS

## Communes de France
mapcom <- st_read("tmp/maps/COMMUNE.shp", quiet=T)

## Aire urbaine
AA_adm<-readRDS(paste0(myrep,"AA_adm.RDS"))
AA_mapcom <- mapcom %>% filter(INSEE_COM %in% AA_adm$CODGEO)
saveRDS(AA_mapcom,paste0(myrep,"AA_mapcom.RDS"))

## EPCI
EPCI_adm<-readRDS(paste0(myrep,"EPCI_adm.RDS"))
EPCI_mapcom <- mapcom %>% filter(INSEE_COM %in% EPCI_adm$CODGEO)
saveRDS(EPCI_mapcom,paste0(myrep,"EPCI_mapcom.RDS"))

## Commune centre
CTR_adm<-readRDS(paste0(myrep,"CTR_adm.RDS"))
CTR_mapcom <- mapcom %>% filter(INSEE_COM %in% CTR_adm$CODGEO)
saveRDS(CTR_mapcom,paste0(myrep,"CTR_mapcom.RDS"))

AA_mapcom<-readRDS(paste0(myrep,"AA_mapcom.RDS"))
EPCI_mapcom<-readRDS(paste0(myrep,"EPCI_mapcom.RDS"))
CTR_mapcom<-readRDS(paste0(myrep,"CTR_mapcom.RDS"))
plot(AA_mapcom$geometry, col="lightyellow", main="zone d'étude")
plot(EPCI_mapcom$geometry,col="orange", add=T)
plot(CTR_mapcom$geometry,col="red", add=T)

## Communes de France
mapiris <- st_read("tmp/maps/CONTOURS-IRIS.shp", quiet=T)


## EPCI
EPCI_adm<-readRDS(paste0(myrep,"EPCI_adm.RDS"))
EPCI_mapiris <- mapiris %>% filter(INSEE_COM %in% EPCI_adm$CODGEO)
saveRDS(EPCI_mapiris,paste0(myrep,"EPCI_mapiris.RDS"))

## Commune centre
CTR_adm<-readRDS(paste0(myrep,"CTR_adm.RDS"))
CTR_mapiris <- mapiris %>% filter(INSEE_COM %in% CTR_adm$CODGEO)
saveRDS(CTR_mapiris,paste0(myrep,"CTR_mapiris.RDS"))

EPCI_mapiris<-readRDS(paste0(myrep,"EPCI_mapiris.RDS"))
CTR_mapiris<-readRDS(paste0(myrep,"CTR_mapiris.RDS"))
plot(EPCI_mapiris$geometry, col="lightyellow", main="EPCI",border="gray50", lwd=0.4)
plot(CTR_mapiris$geometry,col="orange", add=T, border="gray50", lwd=0.4)
plot(EPCI_mapcom$geometry,col=NA, border="black",lwd=1,add=T)

Etape 3 : logements ordinaires en 2019

Nous partirons des fichiers détail de l’INSEE car, à la différence des tableaux prédéfinis, ils permettent virtuellement toutes les formes de croisement d’indicateurs. Ils sont évidemment beaucoup plus volumineux, mais ce sera justement l’occasion pour les étudiants en data mining d’être confrontés à des problèmes d’optimisation et de big data. On trouve leur description détaillée sur le site de l’INSEE

Compte-tenu de la taille des fichiers, nous travaillerons sur une partition de la France en zones

• Zone A : Région Île-de-France (région 11) ;
• Zone B : Régions Centre-Val de Loire (région 24), Bourgogne-Franche-Comté (région 27), Normandie (région 28) et Hauts-de-France (région 32) ;
• Zone C : Régions Grand Est (région 44), Pays de la Loire (région 52) et Bretagne (région 53);
• Zone D : Régions Nouvelle-Aquitaine (région 75) et Occitanie (région 76) ;
• Zone E : Régions Auvergne-Rhônes-Alpes (région 84), Provence-Alpes-Côte d’Azur (région 93), Corse (région 94), Guadeloupe (région 01), Martinique (région 02), Guyane (région 03) et La Réunion (région 04).

# Lit le fichier
men<-fread(ficmen)


# Ajoute la variable INSEE_COM en format caractère en rétablissant le zéro
code <-as.character(men$COMMUNE)
code[nchar(code)==4]<-paste0("0",code[nchar(code)==4])
men$INSEE_COM<-code

## Aire urbaine
AA_adm<-readRDS(paste0(myrep,"AA_adm.RDS"))
AA_men <- men %>% filter(INSEE_COM %in% AA_adm$CODGEO)
saveRDS(AA_men,paste0(myrep,"AA_men.RDS"))

## EPCI
EPCI_adm<-readRDS(paste0(myrep,"EPCI_adm.RDS"))
EPCI_men <- men %>% filter(INSEE_COM %in% EPCI_adm$CODGEO)
saveRDS(EPCI_men,paste0(myrep,"EPCI_men.RDS"))

## Commune centre
CTR_adm<-readRDS(paste0(myrep,"CTR_adm.RDS"))
CTR_men <- men %>% filter(INSEE_COM %in% CTR_adm$CODGEO)
saveRDS(CTR_men,paste0(myrep,"CTR_men.RDS"))

# Lecture du fichier de métadonnées
meta<-fread(ficmen_meta)

# Elimine les modalités IRIS, COMMUNE et TRIRIS
meta <-meta %>% filter(!(COD_VAR %in% c("COMMUNE", "IRIS","TRIRIS")))

# Enregistrement dans le dossier data
saveRDS(object = meta,
        file = paste0(myrep,"men_meta.RDS"))

Le format sf (spatial features)

La cartographie et plus généralement les opérations géométriques sur des données spatiales dans R peuvent facilement être effectuées avec le package sf (spatial features) qui crée des objets uniques rassemblant à la fois

• un tableau de données (l’équivalent du fichier .dbf)
• une géométrie (l’équivalent du fichier .shp)
• une projection (l’équivalent du fichier .prj)

Lorsqu’on récupère des fonds de carte au format shapefile (.shp) ou dans d’autres formats standards comme GeoJson, la première tâche consiste donc à les convertir au formt sf afin de pouvoir les utiliser facilement dans R. L’importation se fait à l’aide de l’instruction st_read en indiquant juste le nom du fichier .shp à charger. Les autres fichiers (.dbf ou .proj) seront lus également et intégrés dans l’objet qui hérite de la double classe data.frame et sf.

Etapes de préparation des données

Dans notre exemple, nous allons suivre les étapes suivantes :

1. Préparer les données statistiques par IRIS dans un data.frame
2. Charger un fonds de carte par IRIS au format sf
3. Effectuer une jointure entre les deux fichiers par le code IRIS
4. Sauvegarder le résultat
5. Agréger les données statistiques et géométriques par commune
6. Sauvegarder le résultat.

Préparer les données statistiques

On importe le fichier des individus et on corrige le code IRIS des communes non divisées en IRIS.

myrep<-"data/Rennes/"
tab_ind<-readRDS(paste0(myrep,"EPCI_men.RDS"))
tab_ind$IRIS[tab_ind$IRIS=="ZZZZZZZZZ"]<-
  paste0(tab_ind$INSEE_COM[tab_ind$IRIS=="ZZZZZZZZZ"],"0000")

head(tab_ind[,1:5],3)

FALSE    COMMUNE   ARM      IRIS ACHL AEMM
FALSE 1:   35001 ZZZZZ 350010101  B12 1976
FALSE 2:   35001 ZZZZZ 350010101  B12 1977
FALSE 3:   35001 ZZZZZ 350010101  B12    0

Selectionner les lignes et colonnes

On décide de ne sélectionner que les ménages de propriétaires installés depuis moins de 3 ans et on retient comme variable le type de logement que l’on recode en deux catégories

tab_ind2 <- tab_ind %>% filter(as.numeric(ANEM) < 3,        # durée d'occupation < 3 ans
                               STOCD =="10",                # propriétaire de son logement
                               TYPL %in% c(1,2)) %>%        #  Maison ou Appartement  
                        mutate(TYPL = as.factor(TYPL)) %>%
                        select(IPONDL, IRIS, TYPL) 
levels(tab_ind2$TYPL)<-c("Maison","Appt")

knitr::kable(head(tab_ind2,5),digits=2)

Agréger les données

On commence par créer un tableau long croisant les deux variables et leur effectif pondéré :

tab_long<- tab_ind2 %>%
           group_by(IRIS,TYPL)%>%
           summarise(NB=sum(IPONDL))


knitr::kable(head(tab_long,5),digits=2)

Pivoter le tableau

Puis on fait “pivoter” le tableau pour l’obtenir en format large :

tab_large <- tab_long %>% pivot_wider(id_cols = IRIS, 
                                      names_from = TYPL,
                                      names_prefix = "T_",
                                      values_from = NB,
                                      values_fill = 0)

knitr::kable(head(tab_large,5),digits=2)

Ajouter de nouvelles variables

On ajoute de nouvelles variables telles que le nombre total de nouveaux propriétaires et la part des propriétaires de maisons parmi eux.

tab<- tab_large %>% mutate(T_TOT = T_Maison+T_Appt,
                           Maison_pct = 100*T_Maison/T_TOT)

knitr::kable(head(tab,5),digits=c(0,0,0,0,1))

Examiner la distribution statistique

On examine l’histogramme donnant distribution statistique du % de nouveux propriétaires en maison individuelle.

p <- ggplot(tab) + aes (x = Maison_pct) +
                   geom_histogram(breaks = c(0,10,20,30,40,50,
                                             60,70,80,90, 100)) +
                   scale_x_continuous("% de Maisons") +
                   scale_y_continuous("Nombre d'IRIS") +
                   ggtitle(label = "Type de logement des nouveaux propriétaires",
                           subtitle = "Source : INSEE, RP 2019 - Métropole de Rennes")
                            
p

Charger les données géométriques

On importe le fichier des iris de l’agglomération qui est de type sf (spatial features)

map_iris <- readRDS(paste0(myrep,"EPCI_mapiris.RDS"))
map_iris<-map_iris[,c(5,6,3,2)]
names(map_iris)<-c("IRIS","NOM_IRIS","COM","NOM_COM","geometry")

class(map_iris)

FALSE [1] "sf"         "data.frame"

knitr::kable(head(as.data.frame(map_iris)[,1:5],2))

Visualisation du fonds iris avec sf

On peut facilement faire un plot de la colonne geometry du fichier sf

plot(map_iris$geometry,col="lightyellow")

Jointure des données IRIS et du fonds de carte

map_iris_tab<-merge(map_iris,tab,
                   by.x="IRIS",by.y="IRIS",
                   all.x=T,all.y=F)

knitr::kable(head(map_iris_tab,3),digits=2)

Cartographie rapide avec sf

Une astuce pour visualiser rapidement une carte choropléthe avec le seul packahe sf :

plot(map_iris_tab[,"Maison_pct"])

Sauvegarde du fichier par IRIS

saveRDS(map_iris_tab,paste0(myrep,"map_iris_typl.RDS"))

st_write(map_iris_tab,paste0(myrep,"map_iris_typl.shp"))

Agrégation statistique + géométriques

Grâce aux nouveaux packages de R (dplyr et sf) il est possible d’agréger simultanément les statistiques et les géométries après les avoir stockés dans un même objet de type “sf”

Du coup, on peut gagner beaucoup de temps dans les traitements et les analyses cartographiques, en particulier si l’on veut tester différents niveaux d’agrégation.

Agrégation des IRIS en communes

L’agrégation est très facile et elle concerne à la fois les variables (de stock) et les geometries

map_com_tab <- map_iris_tab %>% 
  group_by(COM, NOM_COM) %>% 
  summarise(T_Maison=sum(T_Maison,na.rm=T), 
            T_Appt=sum(T_Appt,na.rm=T)) %>%
  st_cast("MULTIPOLYGON")

map_com_tab <- map_com_tab %>%  mutate(T_TOT = T_Maison+T_Appt,
                                  Maison_pct = 100*T_Maison/T_TOT) 

Agrégation des iris en communes

knitr::kable(head(st_drop_geometry(map_com_tab),3),digits=c(0,0,0,0,0,1))

plot(map_com_tab$geometry,col ="lightyellow")

Examiner la distribution statistique

On examine l’histogramme donnant distribution statistique du % de nouveaux propriétaires en maison par commune.

p <- ggplot(map_com_tab) +aes (x = Maison_pct) +
                   geom_histogram(breaks = c(0,10,20,30,40,50,
                                             60,70,80,90, 100)) +
                   scale_x_continuous("% de Maisons") +
                   scale_y_continuous("Nombre de communes") +
                   ggtitle(label = "Type de logement des nouveaux propriétaires",
                           subtitle = "Source : INSEE, RP 2019, Agglomération de Rennes")
p                            

Cartographie rapide

plot(map_com_tab[,"Maison_pct"])

Le package map_sf

Le package mapsf permet de réaliser des cartes statiques de très haute qualité. Il a en effet été mis au point par des cartographes et des géomaticiens professionnels de l’UMS RIATE. Il prend la suite du package cartography dont la maintenance demeurera assuré quelque temps encore mais ne fera plus l’objet de développements futurs. Le package mapsf présente l’avantage d’être totalement compatibvle avec le package sf ce qui n’était pas autant le cas pour le package cartography, plus ancien, et créé pour être compatible avec l’ancien package sp.

On trouvera la documentation du package mapsf à l’adresse suivante :

https://riatelab.github.io/mapsf/index.html

Création d’un template cartographique

Nous allons dans un premier temps apprendre à créer un fonds de carte vierge mais comportant tout l’habillage nécessaire (“template”). Pour cela nous allons charger différentes couches cartographiques correspondant respectivement au département, aux communes et aux iris :

myrep<-"data/Rennes/"
map_iris<-readRDS(paste0(myrep,"map_iris_typl.RDS"))
map_com<-readRDS(paste0(myrep,"map_com_typl.RDS"))

tracé d’un fonds de carte vierge

La fonction mf_map() avec le paramètre type = "base"permet de tracer une carte vide

 mf_map(map_iris, type = "base")

Superposition de couches

On peut toutefois ajouter toute une série de paramètres supplémentaire (col=, border=, lwd=, …) et superposer plusieurs fonds de carte avec le paramètre add = TRUE. L’ajout de la fonction layout permet de rajouter un cadre une légende.

# Trace les Iris avec des paramètres
mf_map(map_iris,  
       type = "base", 
       col = "lightyellow", 
       border="gray50",
       lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com,  
       type = "base", 
       col = NA,
       border="red",
       lwd=0.6,
       add = TRUE)

# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et IRIS de la zone d'étude", 
          credits = "Sources : IGN et INSEE")

Ajout d’un thème

On peut finalement modifier l’ensemble de la carte en lui ajoutant une instruction mf_theme() qui peut reprendre des styles existants ( “default”, “brutal”, “ink”, “dark”, “agolalight”, “candy”, “darkula”, “iceberg”, “green”, “nevermind”, “jsk”, “barcelona”) mais aussi créer ses propres thèmes

#Choix du thème
mf_theme("darkula")
# Trace les Iris avec des paramètres
mf_map(map_iris,  
       type = "base", 
       col = "lightyellow", 
       border="gray50",
       lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com,  
       type = "base", 
       col = NA,
       border="red",
       lwd=0.6,
       add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et IRIS de la zone d'étude", 
          credits = "Sources : IGN et INSEE")

Ajout de texte

On peut ajouter une couche de texte avec la fonction mf_label(). Par exemple, on va ajouter à la carte précédente le nom des communes

mf_theme("agolalight")

# Trace les Iris avec des paramètres
mf_map(map_iris,  
       type = "base", 
       col = "lightyellow", 
       border="gray50",
       lwd=0.3)

# Ajoute les contours des communes
mf_map(map_com,  
       type = "base", 
       col = NA,
       border="red",
       lwd=0.6,
       add = TRUE)

# Ajoute une couche de labels
mf_label(map_com, 
         var="COM",
         cex=0.4,
         col="blue",
         overlap = FALSE)

# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et IRIS de la zone d'étude", 
          credits = "Sources : IGN et INSEE")

Carte de stock

Une carte de stock représente la localisation de quantités que l’on peut aditionner et dont le total a un sens. Par exemple un nombre d’habitants, un nombre de ménages, un nombre d’automobiles. Ce quantités doivent être représentées par des figures (cercles, carrés, …) dont la surface est proportionelle au stock afin que l’oeil du lecteur puisse les aditionner visuellement.

Dans le package mapsf, on réalise ce type de carte à l’aide de la fonction mf_map()en lui donnant le paramètre type="prop".

On va tenter à titre d’exemple de représenter la distribution du nombre de nouveaux propriétaires.

Carte de stock minimale

Les instructions minimales sont les suivantes :

# Trace les contours des communes
mf_map(x= map_iris, 
       type = "base")

# Ajoute le nombre de ménages par IRIS
mf_map(x =map_iris, 
      type ="prop",
      var = "T_TOT",
      add=TRUE)

Mais le résultat est peu satisfaisant car les cercles sont trop grands. Il faut en pratique toujours effectuer un réglage de ceux-ci avec l’instruction inches=

Carte de stock habillée

mf_theme("agolalight")
mf_map(map_iris, type = "base",  
       col = "lightyellow",border="gray80", lwd=0.3)
mf_map(map_com, type = "base", 
       col = NA,border="black",lwd=1,add = TRUE)

mf_map(map_iris, var = "T_TOT",type = "prop",
  inches = 0.05, col = "red",leg_pos = "left",  
  leg_title = "Nombre de ménages", add=TRUE)

mf_layout(title = "Distribution des nouveaux propriétaires", 
          frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte choroplèthe

Une carte choroplèthe ou d’intensité représente un phénomène relatif dont la somme n’a pas de sens. Par exemple, il serait absurde d’aditionner les % de logement HLM des IRIS du Val de Marne. Ces variables d’intensité caractèrisent donc l’état général d’une zone (choros) et elles vont être représentées par une couleur appliquée à toute la surface de la zone, d’où leur nom de cartes choroplèthes.

La fonction du package mapsf adaptée aux variables d’intensité est la fonction mf_map()munie du paramètre type = "choro".

On va prendre l’exemple du % de maisons parmi les nouveaux propriétaires

Carte choroplèthe minimale

Si on ne précise rien, la carte est réalisée à l’aide de la palette par défaut avec un découpage des classes en quantiles (effectifs égaux).

# Carte choroplèthe
mf_map(
  x = map_iris, 
  var = "Maison_pct",
  type = "choro")

Carte choroplèthe habillée

On peut arriver à une carte beaucoup plus satisfaisante en contrôlant l’ensemble des paramètres de couleur et de découpage des classes. Puis en superposant les contours de communes au dessus de la carte des IRIS pour faciliter le repérage.

mybreaks = c(0, 10,20,30,40,50,
             60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5), 
                    pal = c("Greens", "Reds"))
# Carte choroplèthe des iris
mf_map( map_iris, var = "Maison_pct",
        type = "choro",
        breaks = mybreaks,pal = mypal, 
        border=NA,
       col_na = "gray80",leg_title = "% maisons",
       leg_val_rnd = 0)
# Contour des communes et cadre
mf_map(map_com, type = "base", col = NA,
       border="black",lwd=1,add = TRUE)
mf_layout(title = "Les nouveaux propriétaires de maisons", 
          frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte stock + choroplèthe (1)

On peut combiner les deux modes cartographiques par superposition :

mf_theme("agolalight")

# Choisit les classes
mybreaks = c(0,20,40,60,80,100)

# Trace la carte choroplèthe
mf_map(
  x = map_iris, 
  var = "Maison_pct",
  breaks = mybreaks,
 # pal=mypal,
 type = "choro",
  border="white",
  col_na = "gray80",
 lwd=0.3,
 leg_title = "% Maisons", 
 leg_val_rnd = 0,
  
)

# Ajoute les cercles proportionnels

mf_map(
  x =map_iris, 
  var = "T_Maison",
  type = "prop",
  inches = 0.06, 
  col = "red",
  leg_pos = "right",  
  leg_title = "Effectif",
  add=TRUE
)
# Ajoute les contours des communes
mf_map(map_com, 
       type = "base", 
       col = NA,
       border="black",
       lwd=1,
       add = TRUE)

# Ajoute un cadre, un titre et des sources
mf_layout(title = "Les nouveaux propriétaires de maison", 
          frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte stock + choroplèthe (2)

Mais les cercles dissimuent alors les plages de couleur, aussi on peut utiliser le type prop_choro qui place la variable choroplèthe à l’intérieur des cercles

mf_theme("agolalight")
mybreaks = c(0, 10,20,30,40,50,60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5), pal = c("Greens","Reds"))
mf_map(map_iris, type = "base",  
       col = "gray80",border="white", lwd=0.3)
mf_map(map_com, type = "base", 
       col = NA,border="white",lwd=1,add = TRUE)
mf_prop_choro( x = map_iris,  var = c("T_TOT", "Maison_pct"), 
  inches = 0.06, col_na = "grey", pal=mypal,
  breaks = mybreaks, nbreaks = 4, lwd = 0.1,
  leg_pos = c("right", "left"),leg_val_rnd = c(0,0),
  leg_title = c("Nouveaux propriétaires", "% Maison"),
  add = TRUE)
mf_layout(title = "Les nouveaux propriétaires de maisons", 
          frame = TRUE, credits = "Sources : IGN et INSEE")

Introduction

# Commune centre 
codectr<-  "35238"  # Code
namectr <- "Rennes"  # Nom
latctr <-  -1.67789
lonctr <-  48.11204
# Choix du rayon de collecte des dvf (en mètres)
rayon <- 50000
# Dossier de stockage
myrep <- "data/Rennes/"

Données dvf

Les données sur les dvf sont disponibles sur le site publicopendatasof:

https://public.opendatasoft.com

On s’intéresse plus particulièrement aux données dvf géolocalisées accessibles ici :

https://public.opendatasoft.com/explore/dataset/buildingref-france-demande-de-valeurs-foncieres-geolocalisee-millesime/

Sélection des informations

Le fichier comporte 28 millions d’enregistrement ce qui est évidemment beaucoup … On va donc procéder à une sélection en se servant des différents onglets disponibles.

A titre d’exemple, nous allons essayer de télécharger les enregistrements vérifiant les conditions suivantes :

• localisation dans un rayon de 60 km autour d’Amiens
• type de transaction : Vente
• type de bien : Maison ou appartement

Nous constatons qu’il y a 133294 enregistrements vérifiant ces conditions.

Lien de téléchargement

Nous pouvons récupérer les données soit au format .csv pour les analyse statistiques, soit au format .geojson pour les analyses spatiales. Mais dans ce cas on ne peut pas dépasser la taille de 50 000 enregistrement. Il vaut donc mieux télécharger au format .csv puisque ce fichier comporte les latitudes et longitudes ce qui nous permettra de créer un fichier cartographique a posteriori

On peut effectuer le téléchargement depuis le site web ou bien juste enregistrer le lien de téléchargement et effectuer l’opération dans R.

Pour trouver le lien on passe la souris au dessus du lien “télécharger les données au format .csv / seulement les enregistrements sélectionnées” et on effectue un click droit suivi de “copier le lien”

Construction d’une API

On récupère le lien de téléchargement du fichier .csv :

https://public.opendatasoft.com/api/explore/v2.1/catalog/datasets/buildingref-france-demande-de-valeurs-foncieres-geolocalisee-millesime/exports/csv?lang=fr&refine=nature_mutation%3A%22Vente%22&refine=type_local%3A%22Maison%22&refine=type_local%3A%22Appartement%22&facet=facet(name%3D%22nature_mutation%22%2C%20disjunctive%3Dtrue)&facet=facet(name%3D%22type_local%22%2C%20disjunctive%3Dtrue)&where=(distance(%60geo_point%60%2C%20geom%27POINT(2.24771118722856%2049.91161703449722)%27%2C%2060298.74879666418m))&timezone=Europe%2FParis&use_labels=true&delimiter=%3B

A première vue ceci paraît assez obscur mais on réalise assez vite que cette URL reprend l’ensemble des spécifications envoyées pour sélectionner nos enregistrements. On peut souligner en gras les parties utiles :

• &refine=nature_mutation%3A%22Vente%22&
• &refine=type_local%3A%22Maison%22
• &refine=type_local%3A%22Appartement%22
• &where=(distance(%60geo_point%60%2C%20geom%27
• POINT(2.26968%2049.90807)%27%2C%2060037m%20))

Nous allons donc créer un lien qui extrait automatiquement les dvf en modifiant juste la position du point central et le rayon. On prend cette fois-ci nos paramètres de Rennes

myurl<-paste0("https://public.opendatasoft.com/api/explore/v2.1/catalog/datasets/buildingref-france-demande-de-valeurs-foncieres-geolocalisee-millesime/exports/csv?lang=fr&refine=nature_mutation%3A%22Vente%22&refine=type_local%3A%22Maison%22&refine=type_local%3A%22Appartement%22&facet=facet(name%3D%22nature_mutation%22%2C%20disjunctive%3Dtrue)&facet=facet(name%3D%22type_local%22%2C%20disjunctive%3Dtrue)&where=(distance(%60geo_point%60%2C%20geom%27POINT(", latctr,"%20",lonctr,")%27%2C%20",rayon,"m))&timezone=Europe%2FParis&use_labels=true&delimiter=%3B")

Téléchargement

On lance la requête avec l’instruction download.file() et on stocke le résultat dans notre répertoire de travail :

download.file(url = myurl, 
              destfile = paste0(myrep,"dvf.csv"))

essai de l’URL ‘https://public.opendatasoft.com/…’ downloaded 38.4 MB

L’opération s’execute en un peu moins d’une minute dans cet exemple. Mais il peut arriver qu’elle échoue lorsque le serveur est saturé et dans ce cas il faut recommencer.

Introduction

# Commune centre 
codectr<-  "35238"  # Code
namectr <- "Rennes"  # Nom
latctr <-  -1.67789
lonctr <-  48.11204
# Choix du rayon de collecte des dvf (en mètres)
rayon <- 50000
# Dossier de stockage
myrep <- "data/Rennes/"
# Fichier dvf brut
myfic <-"dvf.csv"

Chargement du jeu de données

On adapte le programme de Boris et Florent en utilisant la fonction fread() du package data.table car elle est plus rapide/

DVF<-fread(paste0(myrep,myfic))

Etape 1 : Sélection des mutations de type “Ventes” de “Maison” et “Appartement’

etape1 <- DVF %>% filter(`Nature de la mutation` == "Vente")
etape1bis <- etape1 %>% filter(`Type de local` == "Maison" |
                              `Type de local` == "Appartement")

Etape 2 : Sélection et renommage des variables

etape2 <- etape1bis %>% 
          select(id = `Identifiant de mutation (Etalab)`,
                 disposition = `Numéro de disposition`  , 
                 parcelle= `Identifiant de la parcelle cadastrale`,
                 date = `Date de la mutation`,
                 nature = `Nature de la mutation`, 
                 INSEE_COM = `Code INSEE de la commune`,
                 INSEE_DEP =  `Code INSEE du département` , 
                 type = `Type de local` ,
                 surface = `Surface réelle du bâti` , 
                 piece = `Nombre de pièces principales`, 
                 prix = `Valeur foncière`  , 
                 latitude = `Latitude`, 
                 longitude = `Longitude`)

Etape 3 : Remplacement des cellules vides par des NA et suppression des NA

etape2[etape2==""] <- NA
etape3 <- etape2 %>% na.omit()

Etape 4 : Sélections des mutations simples

unique <- etape3 %>% distinct(id, prix, surface)
nbunique <- unique %>% group_by(id) %>% summarise(nb = n())

etape4 <- nbunique %>% filter(nb==1)

Etape 5 : Jointure attributaire pour récupérer les informations de la mutation

merge <- merge(etape4,etape3, by="id")
etape5 <- merge %>% distinct(id, .keep_all = TRUE) %>% 
select(id, date, type, nature, INSEE_COM, prix,
      surface, piece, latitude, longitude)

etape5$prix <- as.numeric(etape5$prix)
etape5$surface <- as.numeric(etape5$surface)
etape5$piece <- as.numeric(etape5$piece)

quantile(etape5$prix, c(0.01, 0.99))

    1%    99% 
 17000 616320 

hist(etape5$prix, breaks = 50000, xlim = c(15000,1000000))

table(etape5$type)

Appartement      Maison 
      42003       68734 

Maisons <- etape5 %>% filter(type == 'Maison')
Appartement <- etape5 %>% filter (type == 'Appartement')    

par(mfrow=c(1,2))
hist(Maisons$surface, nclass = 500, xlim = c(0,600))
hist(Appartement$surface, nclass = 500, xlim = c(0,300))

Etape 6 : Sélection des bornes de prix et de surface

etape6 <- etape5 %>% 
       filter(between(prix, 15000, 5000000)) %>%
       filter(case_when(type=='Appartement' ~  
                          between(surface, 10, 200)) | 
       case_when(type=='Maison' ~ 
                   between(surface, 10, 350)))

Etape 7 : Création de la variable prix/m2

etape7 <- etape6 %>% mutate(prixm2 = prix/surface)

Etape 8 : Sélection des bornes de prix au m2

quantile(etape7$prixm2,c(0.01,0.99))

       1%       99% 
 387.3057 5825.8048 

hist(etape7$prixm2, breaks = 1000, xlim = c(0,10000))

Filtrer les valeurs aberrantes

etape8 <- etape7 %>% 
          filter(between(prixm2,330,8000))

Etape 9 : Corrections diverses

etape8$date <- as.character(etape8$date)
etape8 <- etape8 %>% 
            mutate(annee = substr(etape8$date, 1, 4))

etape8$prix <- round(etape8$prix)
etape8$prixm2 <- round(etape8$prixm2)

code <-as.character(etape8$INSEE_COM)
code[nchar(code)==4]<-paste0("0",code[nchar(code)==4])
etape8$INSEE_COM<-code
etape8$INSEE_DEP<-substr(etape8$INSEE_COM,1,2)

Etape 10 : Exportation

DVFOK <- etape8 %>% select(id, date, annee = annee, 
                           INSEE_COM, INSEE_DEP,
                           type, prix, surface, prixm2,
                           latitude, longitude)%>%
                     arrange(date)

write.csv(DVFOK, 
          paste0(myrep, "dvf_clean.csv"),
          quote = FALSE)

saveRDS(DVFOK, paste0(myrep, "dvf_clean.RDS"))

Bilan

dim(DVF)[1]

[1] 151749

dim(DVFOK)[1]

[1] 108695

100*(dim(DVF)[1]-dim(DVFOK)[1])/dim(DVF)[1]

[1] 28.37185

Retour sur sf

Nous revenons sur le package sf (spatial features) que nous avons déjà rencontré au moment de la création de cartes thématiques par IRIS ou communes à l’aide du package mapsf.

Ici le package sf va être utilisé pour cartographier non pas des zones mais des localisations ponctuelles. Il pourra être à nouveau couplé avec le logiciel de cartogaphie statique comme mapsf , afin par exemple de placer les localisations des ventes foncières au dessus du fonds de carte des IRIS ou communes.

Mais il pourra aussi servir de base à des cartographies dynamiques permettant de placer les points correspondant aux ventes sur des réseaux de rue et plus généralement sur des “tuiles” cartographiques permettant d’effectur des zoom. On utilisera à cet effet d’autres packages comme leaflet ou sa version simplifiée mapview.

Zone d’étude

On définit les paramètres de la zone d’étude en ajoutant le lien vers les dvf nettoyées

# Commune centre 
codectr<-  "35238"  # Code
namectr <- "Rennes"  # Nom
latctr <-  -1.67789
lonctr <-  48.11204
# Choix du rayon de collecte des dvf (en mètres)
rayon <- 50000
# Dossier de stockage
myrep <- "data/Rennes/"
# Fichier dvf nettoyé
mydvf <-"dvf_clean.RDS"

Choix de communes

On choisit plus précisément de travailler sur les communes de l’agglmomération

map_com <- readRDS(paste0(myrep,"EPCI_mapcom.RDS")) %>%
           select(INSEE_COM, NOM, geometry)
plot(map_com$geometry,
     col="lightyellow")

Sélection des dvf nettoyées

Nous reprenons le fichier nettoyé de localisation des ventes de maison et d’appartement établi au chapitre précédent et nous ne conservons que 9 variables ainsi que les ventes situées dans les communes de la zone choisie

dvf <- readRDS(paste0(myrep,mydvf)) %>%
        select(annee, INSEE_COM,
               type, prix, surface, prixm2,
               latitude, longitude) %>%
        filter(INSEE_COM %in% map_com$INSEE_COM)


kable(head(dvf,3))

Transformation des dvf en fichier sf

On ajoute une colonne geometry en utilisant la fonction st_as_sf()du package sf. On attribue à la projection le code EPSG = 4326 qui correspond à un fichier latitude longitude.

map_dvf <- st_as_sf(dvf, coords = c("longitude","latitude"))
st_crs(map_dvf)<- 4326

kable(head(map_dvf,3))

Test de superposition

On tente de superposer les deux cartes des ventes et des communes.

plot(map_com$geometry, col="lightyellow")
plot(map_dvf$geometry, col="red", cex=0.2, pch=20,add=T)

plot(map_dvf$geometry, col="red", cex=0.2, pch=20)
plot(map_com$geometry, col="lightyellow", add=T)

Vérification de la projection

Nous savons que les coordonnées du fichier DVF sont en latitude longitude (EPSG=4326). Mais la projection de notre fonds des communes est en Lambert-93 (EPSG = 2154)

st_crs(map_com)

Coordinate Reference System:
  User input: RGF93 v1 / Lambert-93 
  wkt:
PROJCRS["RGF93 v1 / Lambert-93",
    BASEGEOGCRS["RGF93 v1",
        DATUM["Reseau Geodesique Francais 1993 v1",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4171]],
    CONVERSION["Lambert-93",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",46.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",3,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",49,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",44,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",700000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",6600000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["France - onshore and offshore, mainland and Corsica."],
        BBOX[41.15,-9.86,51.56,10.38]],
    ID["EPSG",2154]]

A priori il s’agit bien de la même de sorte que les coordonnées X,Y devraient bien se superposer sur le fonds IRIS

Ajustement des deux projections

On transforme donc le fichier des DVF en Lambert-93 pour obtenir une adéquation avec le fichier des communes.

map_dvf <- st_transform(map_dvf,2154)

Et cette fois-ci on réussit à superposer les deux fonds de carte

plot(map_com$geometry, col="lightyellow")
plot(map_dvf$geometry, col="red", cex=0.2, pch=20,add=T)

Cartographie ponctuelle avec mapsf

On peut utiliser les fonctions de mapsf vues précédemment. Par exemple repérer les ventes de maison et d’appartement :

mf_map(map_com, type="base",
       col="lightyellow")
mf_map(map_dvf, type = "typo",
       var= "type",add =T,
       border=NA,
       cex=0.3,
       leg_title = "Bien vendu")
mf_label(map_com,
         var="NOM",cex = 0.4,
         overlap = F,
         col="black",
         halo=T)
mf_layout(title = "Ventes de maisons et d'appartement 2013-2020",
          credits = "Sources : IGN, INSEE, DVF",
          scale = T)

Cartographie ponctuelle avec mapview

Le package mapview pemet de créer facilement des cartes interactives ù l’on pourra visualiser les ventes sur un fonds de carte de son choix (OSRM, ESRI, …). On peut l’utiliser par exemple pour voirune commune précise

library(mapview)
map_dvf_sel<-map_dvf %>% filter(INSEE_COM == "35196")
mapview(map_dvf_sel)

On peut améliorer sensiblement la visualisation à l’aide de paramètres présentés sur le site web du package https://r-spatial.github.io/mapview/

library(mapview)
# Choix de la zone d'étude
map_com_sel<-map_com %>% filter(INSEE_COM == "35196")
map_dvf_sel<-map_dvf %>% filter(INSEE_COM == "35196")

# Choix des tuiles
mapviewOptions(basemaps = c("OpenStreetMap" ,
                           "Esri.WorldImagery"))
# Première couche
map1 <- mapview(map_com_sel,
                col.regions="lightyellow",
                alpha.regions = 0.4)

# Deuxième couche
map2<-mapview(map_dvf_sel,
              zcol = "type",
              cex = 4,
              col.regions=c("blue","red"),
             alpha.regions = 0.4)
# Assemblage
map1+map2

Agrégation par communes

On peut évidemment agréger les informations sur les valeurs foncières par communes. Par exemple compter le nombre de maisons vendues et leur prix moyen au m2

dvf_by_com <- dvf %>% 
               group_by(INSEE_COM) %>%
               filter(type == "Maison") %>%
               summarise(nb = n(),
                         med_prixm2 = median(prixm2))

kable(head(dvf_by_com,3))

On peut désormais effectuer la jointure entre les données agrégées par adresse et le fichier sf de localisation des adresses :

map_com_dvf <- inner_join(map_com,dvf_by_com) %>% st_as_sf()

On peut désormais utiliser les méthodes de cartographie déjà vues avec mapsf :

mf_theme("agolalight")
mybreaks = c(1000, 1500, 2000, 2500, 3000, 3500,
             4000, 4500, 5000, 5500, 6000)
mypal=brewer.pal(n = 10,name = "Spectral")
mf_map(map_com_dvf, type = "choro", 
       var = "med_prixm2", 
       breaks=mybreaks, pal = mypal,
       border="gray", lwd=0.3,
       leg_title = "prix en €/m2", leg_val_rnd = 0,
       leg_pos = "left")

mf_map( map_com_dvf,  type = "prop",
        var = "nb", inches = 0.12, 
        col="black",border="white",
        leg_title = "nb. de ventes",
        leg_pos = "topleft")

mf_layout(title = "Les ventes de maison de 2013 à 2020",
        frame = TRUE, credits = "Sources : IGN et DVF", 
        arrow = F )

Cartographie par IRIS

Supposons maintenant que l’on souhaite cartographier les prix des maisons par IRIS à l’intérieur de la commune de Rennes. Comment faire dans la mesure où l’IRIS n’est pas mentionné dans le fichier DVF ? Réponse : en superposant les deux fonds de carte et en les intersectant.

map_iris_ctr<-readRDS(paste0(myrep,"CTR_mapiris.RDS")) %>%
               select("NOM_COM","INSEE_COM","CODE_IRIS","NOM_IRIS","geometry")
map_dvf_ctr <- map_dvf %>% filter(INSEE_COM %in% map_iris_ctr$INSEE_COM)
plot(map_iris_ctr$geometry, col="lightyellow")
plot(map_dvf_ctr$geometry,col="red",pch=21,cex=0.1,add=T)

inter<-st_intersection(map_dvf_ctr,map_iris_ctr)
kable(head(inter,4))

map_iris_dvf <- inter %>% 
               st_drop_geometry() %>%
               group_by(CODE_IRIS, NOM_IRIS) %>%
               filter(type == "Maison") %>%
               summarise(nb = n(),
                         med_prixm2 = median(prixm2)) %>%
              filter(nb >=5) %>% 
              right_join(map_iris_ctr) %>% 
              st_as_sf()
               
             
kable(head(map_iris_dvf,3))

mf_theme("agolalight")
mybreaks = c(1000, 1500, 2000, 2500, 3000, 3500,
             4000, 4500, 5000, 5500, 6000)
mypal=brewer.pal(n = 10,name = "Spectral")
mf_map(map_iris_dvf, type = "choro", 
       var = "med_prixm2", 
       breaks=mybreaks, pal = mypal,
       border="gray", lwd=0.3,
       leg_title = "prix en €/m2", leg_val_rnd = 0,
       leg_pos = "left")

mf_map(map_iris_dvf,  type = "prop",
        var = "nb", inches = 0.08, 
        col="black",border="white",
        leg_title = "nb. de ventes",
        leg_pos = "topleft")

mf_layout(title = "Les ventes de maison de 2013 à 2020",
        frame = TRUE, credits = "Sources : IGN et DVF", 
        arrow = F )

Définir le sujet

Soit le sujet : Qui sont les nouveaux propriétaires ? Quel est leur niveau de qualification ?

Formuler des hypothèses

Justes ou fausses, les hypothèses permettent de cadrer l’analyse. Elles peuvent faire intervenir plusieurs variables

Charger les données

On charge ici les données ménages de l’aire urbaine.

tab_ind<-readRDS(paste0(myrep,"AA_men.RDS"))
head(tab_ind[,1:5],2)

FALSE    COMMUNE   ARM      IRIS ACHL AEMM
FALSE 1:   22036 ZZZZZ ZZZZZZZZZ  A12 1973
FALSE 2:   22036 ZZZZZ ZZZZZZZZZ  B11 1992

Préparer l’analyse

• Soit la relation entre statut d’occupation (Y) et diplôme le plus élevé du chef de ménage (X). Il s’agit de deux variables catégorielles (= qualitatives) que l’on va typiquement mettre en relation à l’aide d’un tableau de contingence et d’un test du chi-2. L’analyse statistique est simple sous R mais il faut tenir compte de trois difficultés

• Le choix de la population de référence est important. Ici on va sélectionner les ménages installés depuis moins de 3 ans

• la sélection ou le regroupement des diplômes est également important car cela va influer sur les résultats du test.

• la pondération des individus doit également être prise en compte puisque le recensement est basé sur un sondage dans les zones urbaines

Sélection des individus et des variables

tab_sel<- tab_ind %>% 
  filter(as.numeric(ANEM) < 3) %>%
  select(INSEE_COM, ANEM, DIPLM,STOCD, IPONDL) 


knitr::kable(head(tab_sel,4))

Recodage des modalités

On cherche le code des modalités dans le fichier des métadonnées

meta<-readRDS(paste0(myrep,"men_meta.RDS"))
metasel <- meta %>% 
  filter(COD_VAR %in% c("STOCD", "DIPLM"))
kable(metasel[,c(1,3,4)])

On recode les modalités des deux variables en regroupant certaines.

tab_sel$STOCD<-as.factor(tab_sel$STOCD)
levels(tab_sel$STOCD)<-c("Proprétaire","Locataire",
                         "Locataire","Locataire",NA)
tab_sel$DIPLM<-as.factor(tab_sel$DIPLM)
levels(tab_sel$DIPLM) <- c("Aucun","Aucun","Aucun",
                           "BEP","BEP","BEP", "BAC","BAC",
                         "BAC+123","BAC+123","> BAC+3","> BAC+3",NA)

                        
knitr::kable(head(tab_sel,3))

Tableau de contingence (FAUX)

Le plus simple semble être l’instruction table()

tab_cont<-table(tab_sel$STOCD,tab_sel$DIPLM)

knitr::kable(addmargins(tab_cont))

Tableau de contingence (JUSTE)

On pondère avec wtd.table() du package questionr.

library(questionr)
tab_cont_wtd<-wtd.table(tab_sel$STOCD,tab_sel$DIPLM,
                        weights = tab_sel$IPONDL)

knitr::kable(round(addmargins(tab_cont_wtd),0))

Comparaison des résultats

• Tableau non pondéré … légèrement faux !

• Tableau pondéré … juste !

Visualisation du tableau de contingence

On choisit l’orientation du tableau et on l’affiche avec plot()

mytable<-wtd.table(tab_sel$DIPLM,tab_sel$STOCD,weights = tab_sel$IPONDL)
plot(mytable)

Tant qu’à faire, on améliore la figure avec des paramètres supplémentaires :

plot(mytable, 
     main = "Propriété et diplôme dans l'aire urbaine", 
     sub = "Source : INSEE - RP 2019 - Ménages installés depuis - de 3 ans",  
     col=c("lightyellow","lightgreen"))

Test du Chi-deux

Ce test se réalise facilement sur le tableau de contingence avec l’instruction chisq.test() :

mytest<-chisq.test(mytable)
mytest

    Pearson's Chi-squared test

data:  mytable
X-squared = 1731.7, df = 4, p-value < 2.2e-16

• Interprétation : La relation est très significative (p<0.001)

Prolongements

On peut introduire toute une série de variantes à partir de ce premier programme. Par exemple, si on retire la commune-centre, les résultats sont assez différents :

tab_sel2<- tab_sel %>% 
  filter(INSEE_COM != codectr)

mytable<-wtd.table(tab_sel2$DIPLM,tab_sel2$STOCD,weights = tab_sel2$IPONDL)

plot(mytable, 
     main = "Propriété et diplôme hors commune-centre", 
     sub = "Source : INSEE - RP 2019 - Ménages installés depuis - de 3 ans",  
     col=c("lightyellow","lightgreen"))

chisq.test(mytable)

    Pearson's Chi-squared test

data:  mytable
X-squared = 733.11, df = 4, p-value < 2.2e-16

A vous de deviner pourquoi :-)

Formuler l’hypothèse

Soit l’ensemble des ventes effectuées dans un rayon de 50 km autour d’une agglomération. Existe-t-il une relation entre (Y) le prix de vente mesuré en € par m2 de surface habitable et (X) la distance au point central de l’agglomération ?

N.B. Nous retirons la commune centre de l’analyse car on suppose que son comportement interne obéit à une logique différente

Définir l’espace d’étude

Nous allons retenir la distance à vol d’oiseau entre chaque point ou a eu lieu une vente et le centre de l’agglomération que l’on fixera en fonction de notre connaissance du terrain.

# Commune centre 
codectr<-  "35238"  # Code
namectr <- "Rennes"  # Nom

# Choix du rayon de collecte des dvf (en mètres)
rayon <- 50000

# Dossier de stockage
myrep <- "data/Rennes/"

# Fichier dvf nettoyé
mydvf <-"dvf_clean.RDS"

Calculer les distances au centre

On doit ajouter à notre fichier dvf la distance au point central retenu.

## Fichier des dvf
dvf<-readRDS(paste0(myrep,mydvf))

# Aire urbaine
map_com<-readRDS(paste0(myrep,"AA_mapcom.RDS"))

# Selection des dvf
dvf<-dvf %>% filter(INSEE_COM %in% map_com$INSEE_COM) %>%
             filter(INSEE_COM != codectr)

# Transformation cartographique des dvf
map_dvf <- st_as_sf(dvf,coords =c("longitude","latitude")) 
st_crs(map_dvf)<-4326
map_dvf<-st_transform(map_dvf,2154)

## Commune centre => centroïde
map_ctr <-readRDS(paste0(myrep,"CTR_mapcom.RDS"))
map_ctr <-st_centroid(map_ctr)

## Calcul de la distance
map_dvf$dist<-as.numeric(st_distance(map_dvf,map_ctr))/1000

## Visualisation rapide
plot(map_dvf[,"dist"])

Modèle linéaire

On teste l’existence d’un modèle linéaire en fixant les paramètres suivants :

Y <-map_dvf$prixm2
X <-map_dvf$dist
nomY <- "Prix au m2"
nomX <- "Distance en km"
model<- "Modèle Linéaire : Y = aX+b"

mod1 <- lm(Y~X)
summary(mod1)

Call:
lm(formula = Y ~ X)

Residuals:
    Min      1Q  Median      3Q     Max 
-2309.1  -368.5   -47.0   314.5  5856.8 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2837.7774     5.9202   479.3   <2e-16 ***
X            -45.1457     0.3075  -146.8   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 633.1 on 50618 degrees of freedom
Multiple R-squared:  0.2986,    Adjusted R-squared:  0.2986 
F-statistic: 2.155e+04 on 1 and 50618 DF,  p-value: < 2.2e-16

plot(X,Y,
     xlab = nomX,
     ylab = nomY,
     main = model,
     pch=20,
     col="red",
     cex=0.1,
     sub = " Source : DVF 2013-2020"
     )

abline(mod1, col="blue", lwd=2)

Modèle exponentiel

On teste l’existence d’un modèle exponetielle en fixant les paramètres suivants :

Y <-log(map_dvf$prixm2)
X <-map_dvf$dist
nomY <- "log(Prix au m2)"
nomX <- "Distance en km"
model<- "Modèle Exponentiel : log(Y) = aX+b"

mod2 <- lm(Y~X)
summary(mod2)

Call:
lm(formula = Y ~ X)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.95482 -0.14476  0.02937  0.19630  1.68254 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.9791861  0.0031466  2535.8   <2e-16 ***
X           -0.0245287  0.0001635  -150.1   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3365 on 50618 degrees of freedom
Multiple R-squared:  0.3079,    Adjusted R-squared:  0.3079 
F-statistic: 2.252e+04 on 1 and 50618 DF,  p-value: < 2.2e-16

plot(X,Y,
     xlab = nomX,
     ylab = nomY,
     main = model,
     pch=20,
     col="red",
     cex=0.1,
     sub = " Source : DVF 2013-2020"
     )

abline(mod2, col="blue", lwd=2)

Modèle logarithmique

On teste l’existence d’un modèle logarithmique en fixant les paramètres suivants :

Y <-map_dvf$prixm2
X <-log(map_dvf$dist)
nomY <- "Prix au m2"
nomX <- "log(Distance en km)"
model<- "Modèle Logarithmique : Y = a.log(X)+b"

mod3 <- lm(Y~X)
summary(mod3)

Call:
lm(formula = Y ~ X)

Residuals:
    Min      1Q  Median      3Q     Max 
-2820.5  -364.0   -31.8   320.7  5900.1 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3801.054     12.449   305.3   <2e-16 ***
X           -649.871      4.559  -142.5   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 638.5 on 50618 degrees of freedom
Multiple R-squared:  0.2864,    Adjusted R-squared:  0.2864 
F-statistic: 2.032e+04 on 1 and 50618 DF,  p-value: < 2.2e-16

plot(X,Y,
     xlab = nomX,
     ylab = nomY,
     main = model,
     pch=20,
     col="red",
     cex=0.1,
     sub = " Source : DVF 2013-2020"
     )

abline(mod3, col="blue", lwd=2)

Modèle puissance

On teste l’existence d’un modèle puissance en fixant les paramètres suivants :

Y <-log(map_dvf$prixm2)
X <-log(map_dvf$dist)
nomY <- "log(Prix au m2)"
nomX <- "log(Distance en km)"
model<- "Modèle Logarithmique : Y = a.log(X)+b"

mod4 <- lm(Y~X)
summary(mod4)

Call:
lm(formula = Y ~ X)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.13041 -0.14325  0.04097  0.20474  1.55332 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.467618   0.006720  1260.1   <2e-16 ***
X           -0.339949   0.002461  -138.1   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3447 on 50618 degrees of freedom
Multiple R-squared:  0.2738,    Adjusted R-squared:  0.2738 
F-statistic: 1.908e+04 on 1 and 50618 DF,  p-value: < 2.2e-16

plot(X,Y,
     xlab = nomX,
     ylab = nomY,
     main = model,
     pch=20,
     col="red",
     cex=0.1,
     sub = " Source : DVF 2013-2020"
     )

abline(mod4, col="blue", lwd=2)

Choix du meilleur modèle

On décide ici de retenir le modèle 2 (exponentiel) et de stocker les valeurs estimées et résiduelles des prix au m2.

mod<-mod3
map_dvf$prixm2_estim <- mod$fitted.values
map_dvf$prixm2_resid <- mod$residuals

Agrégation par commune

On agrège les valeurs estimées et résiduelles par commune

mod_com <- st_drop_geometry(map_dvf) %>%
              group_by(INSEE_COM) %>%
              summarise(nb=n(),
                        prixm2 = mean(prixm2),
                        prixm2_estim = mean(prixm2_estim),
                        prixm2_resid = mean(prixm2_resid)) %>%
              arrange(prixm2_resid)

On en déduit les communes où le prix est plus bas que ce que laisserait prévoir la distance au centre :

kable(head(mod_com), digits=0)

Et celles où le prix est plus élevé que ce que laisserait prévoir la distance au centre

kable(tail(mod_com),digits=0)

Cartographie des résultats

En se limitant à l’aire urbaine, on va visualiser les résultats du modèle. On commence par effectuer la jointure entre les résultats du modèle de régression et le fonds de carte des communes :

map_com<-map_com %>%
         left_join(mod_com)

On peut alors cartographier les prix observés et théoriques :

mybreaks <-c(1000,1500, 1750, 2000, 2250, 2500, 2750,  3000, 6000)
mypal <- brewer.pal(n = 8, name="YlOrRd")

par(mfrow = c(1,2))
mf_map(map_com, 
       type = "choro",
       var="prixm2_estim",
       breaks=mybreaks,
       pal=mypal,
       leg_val_rnd = 0,
       col_na = "gray80",
       leg_pos = "topleft")
mf_layout("Prix théoriques", 
          arrow=F, frame=T,
          credits = "Source : DVF, 2013-2020")


mf_map(map_com, 
       type = "choro",
       var="prixm2",
       breaks=mybreaks,
       pal=mypal,
       leg_val_rnd = 0,
        col_na = "gray80",
       leg_pos = "topleft")
mf_layout("Prix observés",
          arrow=F, frame=T,
          credits = "Source : DVF, 2013-2020")

Et en déduire la carte des résidus :

mybreaks <-c(-700, -500, -300,-100,100,300,500,700)
mypal <- brewer.pal(n = 7, name="RdBu")
par(mfrow = c(1,1))
mf_map(map_com, 
       type = "choro",
       var="prixm2_resid",
       breaks=mybreaks,
       pal=mypal,
       col_na = "gray80",
       leg_val_rnd = 0,
       leg_pos = "topleft")
mf_layout("Prix résiduels", 
          arrow=F, frame=T,
          credits = "Source : DVF, 2013-2020")